---
title: "R10.Cadastre_Year"
output: pdf_document
---

# Load Packages
```{r}
library(tidyverse)
library(sandwich)
library(stargazer)
```

# Load Data
```{r}
load("Data/data_final.Rdata")
load("Data/cadastre_points.Rdata")
```

# Run Models (App. Table 3cii)
```{r}

muni_comp <- data_final %>%
  group_by(MFULL) %>% 
  summarise(CA1_CAR1U_MeanF = weighted.mean(CA1_CAR1U_MeanF, Area_km2)) %>% # estimate weighted average of points index by parish size
  left_join(cadastre_points) # merge year of first cadastre

# # Histogram to identify potential outliers
# muni_comp %>% ggplot(aes(x=First, fill = Concelho == "CUBA")) + geom_histogram()

m.1 <- lm(CA1_CAR1U_MeanF ~ First, data = muni_comp)
m.2 <- lm(CA1_CAR1U_MeanF ~ First, data = muni_comp %>% filter(Concelho != "CUBA"))

# Robust standard errors
s.1 <- sqrt(diag(vcovHC(m.1, type = "HC1")))
s.2 <- sqrt(diag(vcovHC(m.2, type = "HC1")))

stargazer(m.1, m.2, 
          type = "text",
          se = list(s.1, s.2), 
          covariate.labels = c("Cadastral Year"), 
          dep.var.labels = "Average Points Index")

```

